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a^ ■ 



The two-dimensional Hubbard model is analyzed in the framework of the two-pole expansion. It 
is demonstrated that several theoretical approaches, when considered at their lowest level, are all 
equivalent and share the property of satisfying the conservation of the first four spectral momenta. 
It emerges that the various methods differ only in the way of fixing the internal parameters and that 
■ it exists a unique way to preserve simultaneously the Pauli principle and the particle-hole symmetry. 

I A comprehensive comparison with respect to some general symmetry properties and the data from 

. quantum Monte Carlo analysis shows the relevance of imposing the Pauli principle. 
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^ , I. INTRODUCTION. 

The Hubbard model [1] plays a key role [2] in the theoretical analysis of highly correlated electron systems and in the 
last years many works have been devoted to its study. The relevance of the model is that contains as a main ingredient 
''O ' a competition between the band energy, which tends to distribute the electrons throughout the entire crystal, and the 
Ch ' repulsive electron-electron interaction which tends to localize the electrons on the lattice sites. Although the model 
^ is oversimplified and in spite of a very intensive study, both analytical and numerical, exact solutions exist only in 
I I one and infinite dimensions. For finite dimensions greater than one the study is far from being completed and more 
\ analysis is required. Many techniques and various approximation schemes have been formulated; among them we 
O^l > recall the Hubbard I approximation [1], the slave boson method [3-5], the non-crossing approximation [6-8], the doo 
^ method [9-12], the projection operator method [13- 18], the method of equation of motion [19-21], the spectral density 
approach [22-24], the composite operator method [25-31]. 
' When the lowest level of approximation is considered, many methods [1, 13-31], though based on different schemes, 
QQ , give an equivalent framework of calculation, where the spectral function of the single-particle propagator is approxi- 
' mated by a two-pole expansion. A common feature of all these methods is that a self-consistent procedure is necessary 
C , in order to calculate some parameters which appear in the expression of the spectral function. As a consequence, 
P^ . although many methods seems to be equivalent, remarkable differences arise, according to the different adopted 
procedures. 

The purpose of this work is to analyze various methods and compare the different predictions with respect to some 
general symmetry properties of the model and with respect to the numerical data obtained by quantum Monte Carlo 
(qMC) techniques. 
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O ■ II. THE MODEL 

■ The Hubbard model is defined by 

rS 

^. H = Y,t,,cH^MJ) + uY,n^{^)n^{^)-^iY.cH^)c{^) (2.1) 

ij i i 

The first term is a kinetic term that describes the motion of the electrons among the sites of the Bravais lattice, 
defined by the vector set {Ri}. For a two-dimensional squared lattice and by restricting the analysis to first nearest 
neighbors, the hopping matrix tij has the form 
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t,j = -Ata,j = -4*1 Yl e'''-^^'-'''^a{k) (2.2) 

k 

a(k) = i [cos(fca;a) + cos(A;j,a)] (2.3) 

a being the lattice constant. In addition to the band term, the model contains an interaction term which approximates 
the interaction among the electrons. In the simplest form of the Hubbard model, the interaction is between electrons 
of opposite spin on the same lattice site; the strength of the interaction is described by the parameter U. Various 
generalizations of the model take into account hopping between second nearest neighbors and intersite interactions. 
c{i), c^{i) are annihilation and creation operators for c-electrons at site i in the spinor notation 



= (c| c\) (2.4) 



ncr{i) = c^^{i)ccr{i) is the number operator of electrons with spin a — (j, i) at the «-th site, /z is the chemical potential 
and is introduced in order to control the band filling n = (c^(i)c(i)). By varying the ratio U/t, the particle number 
n and the temperature T, it is believed that the Hubbard model is capable to describe many properties of strongly 
correlated fermion systems [2]. 

The field c{i) satisfies the equation of motion 

*4c(«) = -Mc(i) + ^ti,c(j) + Un{i)c{i) (2.5) 



where the composite field 



r]{i) = n{i)c{i) (2.6) 
appears. By means of the Hamiltonian (2.1) one can derive the equation of motion for the new field 

ig-^v{i) = -(m - U)rj{i) + (2.7) 



where a higher order composite field appears 

'1 



7r(i) = ^ ai 
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^a''n^ii)c{j) + cii)c\j)cii) 



(2.8) 



In Eq. (2.8) we introduced the charge (/U = 0) and spin (/i = 1, 2, 3) density operator 

n^{i) = c\i)a^,c{i) (2.9) 

with the notation 

a^ = {l,a) a^' = {-l,a) (2.10) 

CTfe being the Pauli matrices. 

In order to close the infinite hierarchy of equations of motion some truncation is necessary. One procedure, common 
to various methods, is to choose a basis of operators {^'(i)} and linearize the equation of motion as 

=Yeij^ii) (2.11) 

j 

where the eigenvalue or energy matrix is self-consistently calculated by means of the equation 

e,j = {mi),H],^\j)U)/{[^{i),^\j)U) (2.12) 

The symbol [. . .]± denotes equal-time anticommutator or commutator, in dependence of the statistics of the set {^'(i)}. 
The rank of the energy matrix is equal to the number of components of \E'(«). It should be noted that Sij may extend 
over several lattice points, according to the choice of the basis {'J'(i)}. Different choices for the basis {'S'(i)} can be 
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Refs. 1,20,21 


* = 


V '7/ 


Ref. 18 


* = 


\r,-{n/2)c) 


Refs. 17,26 


* = 


h 



TABLE I. 



taken. The mean field approximation corresponds to the choice ^ = c. In Table I we summarize some choices used 
in various works. 

The composite operator ^{i) is defined as 



= c(i) - r]{i) = [1 - n{i)\c{i) 



(2.13) 



It is easy to show by direct calculations of the energy matrices that all basis are equivalent and they all give place to 
the same linearized equations of motion. By considering a paramagnetic ground state, the thermal retarded Green's 
function 



S^^{i,j) = {R[^{i)^\j)]) 
after the Fourier transformation, has the following expression 

2 (i) X 

^ uj - Ei{k) + 17] 

The energy spectra Ei (k) are given by 

Ei(k) = R{k) + Q(k) £;2(k)=i?(k)-Q(k) 

where 

E(k) = i [-2/. + U]-^^ [m(k) + 8to(k)Ji72] 



(2.14) 

(2.15) 

(2.16) 
(2.17) 



and the following notation has been used 



/i = 1 - n/2 h = n/2 

m(k) =4t[A + a(k)(p - h)] 



(2.18) 



(2.19) 



5(k) = -t/+i-^m(k) 

The parameters A and p describe a constant shift of the bands and a bandwidth renormalization, respectively. They 
are static intersite correlation functions defined as 



A^(r(i)€^(i))-(ry"(OryHO> 



(2.20) 
(2.21) 



The notation c"(i) stands to indicate the field c on the first neighbor sites: 
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(2.22) 



The explicit expressions of the spectral functions a^\, (k) are given by [26] 



1 + 



2Q(k) 

m(k) 
2Q(k) 



h 



2Q(k) 



2Q(k) 



(2.23) 



4V(k) = 



1 



5(k) 



2Q(k) 



^(2) 
'7') 



(k) 



2 



1 + 



g(k) 
2Q(k) 



The spectral functions relative to other basis can be obtained from (2.23) by recalling that c{i) = + ri{i). 
Straightforward calculations, reported in Appendix B, show that the first four momenta 

/ + 00 
dw^(k,a;)a;" 
-oo 

of the electron spectral density 

1 ^ 
A(k, w) = — /m5ec(k,a;) = Va«(k)^ [a; - ^^(k)] 

^ i=i 

satisfies the exact relations 

M(")(k) =i?^.T.( 



dt 



n—p 



oil). 



d \^ 

I ct(j)[ > , 0<p<n 

J E.T. 



dt' 



(2.24) 



(2.25) 



(2.26) 



where F.T. denotes the Fourier transform. Then, there is a complete equivalence between the scheme of calculation 
traced above and the spectral density approach [23], when the two-pole approximation is considered. 

Summarizing, when the lowest order is considered all various methods [1, 13-31] arc equivalent and they all corre- 
spond to a two-pole approximation with the same expression for the spectral function. The correspondence with the 
notation of Refs. 20 and 18 is given by 



Ji/2iy (k) = hh [Wo + Wi (k)] = -hhii - 4tA + 4to(k) [/| - p] 



(2.27) 



However, the scheme is not complete unless the self-consistent procedure is defined. The spectral functions contain 
three unknown parameters, /i, A, and p very diS^erent results are obtained according to the methods used to calculate 
them. This will be discussed in the next Sections. 



III. SELF-CONSISTENT EQUATIONS. 

In Section 2 the energy matrix and therefore the properties of the basic fields {^} have been fixed by means of Eq. 
(2.12). However, as we have seen, this constrain does not fix in a unique way the dynamics; the parameters /U, A, and 
p remain to be determined and more conditions are required. 

One constrain common to all methods is given by the requirement that the particle number is an external parameter, 
fixed by the boundary conditions. Then, the following self-consistent equation must be satisfied: 

n = 2 [1 - Ccc] (3.1) 

where the correlation function Ccc is defined in (A. 8). We need two more conditions. In Hubbard I approximation [1] 
a simple factorization procedure is used for the two-particle Green's function. Then, from the definitions (2.20) and 
(2.21) one obtains the expressions 
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A = (3.2) 

P = Y (3-3) 

However, Eq. (3.2) is not consistent with the present scheme of calculation. The parameter A is directly connected 
to the single-particle Green's function, and from the definition (2.20) one can immediately derive the self-consistent 
equation 

A = C^^ - C^^ (3.4) 

where the time-independent correlation fimctions are defined by Eq. (A. 6). The parameter p plays an important 
role since it is related to neighboring correlations of the charge, spin and pair. An accurate determination of this 
parameter is very important; as it will be shown in the next Sections, different methods of computing p give rise to 
very different physical solutions. 

In the original work by Roth [20] and in subsequent works [18, 21] the parameter p has been calculated by means 
of the equation of motion in the linearized form (2.11). This procedure leads to the following relation 

where 

/?= ^"^7" ^ (3.6) 

-'l-'2 

and the correlation functions C"^, C"^ are given by Eq. (A. 8). 

In the Composite Operator Method we adopt a different procedure to calculate the parameter p. This quantity is 
not expressed in terms of the single-particle propagator, and there is some freedom in its determination. In COM 
we take advantage of this freedom and we fix the parameter p in such a way that the Hilbert space has the right 
properties to conserve the relations among matrix elements imposed by symmetry laws. 

There is a wide agreement that the unusual and somehow unexpected properties observed in the new class of 
materials are due to the presence of a high correlation among the electrons. If the electron interaction is the key 
to understand these materials, then it is very crucial that the symmetry required by the Pauli principle is correctly 
treated. A convenient way to take care of the Pauli principle is to operate in the representation of second quantization 
where the Pauli principle manifests through the algebra. However, algebra is only one ingredient; physical quantities 
are expressed in terms of expectation values of operators and a suitable choice of the Hilbert space must be made. 
Physical laws, expressed as algebraic relations among the observables, manifest at the level of observation as relations 
among matrix elements. 

In a physics dominated by a high correlation among the electrons, we feel that greater attention should be dedicated 
to the conservation of the Pauli principle. It is well known that in most of the approximation schemes this symmetry 
is violated [19], and some alternative schemes should be considered. The Pauli principle requires that 

^.{iWaii.) = c,(i) [n.^ii) - nl^ii)] 4(z) = (3.7) 

At level of matrix elements, this condition requires that 

C^r, = mW{i))=Q (3.8) 

By recalling the expression (A. 5) for the correlation function and by means of Eqs. (A. 7) and (2.19), we easily 
obtain the following self-consistent equation for the parameter p 

(3.) 

where the quantities Fq and Fx are defined by (A. 7). 

Summarizing, the adopted self-consistent procedures are based on the use of Eqs. (3.1), (3.4), (3.5) in the Roth's 
method, and of Eqs. (3.1), (3.4), (3.9) in the COM. It should be noted that the self- consistent equations are all 
coupled, so that a different choice for the third equation will have influence also on the first two equations. In 
particular, when the Pauli condition (3.8) is not satisfied, there is an ambiguity in writing the first-self consistent 
equation (3.1). In the next section we shall discuss and compare the different results which arise, due to the use of 
different self-consistent schemes. 
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APPENDIX A: EQUAL-TIME CORRELATION FUNCTIONS 



By standard methods equal-time correlation functions can be calculated by means of the retarded Green's function. 
Indeed, one has 



(27r)2 



ttJ 



(A.l) 



where 5ii,*(k, lu) is the Fourier transform of the Green's function, fpi^) is the Fermi distribution function and Q.b is 
the first Brillouin zone. By means of Eq. (2.15) we have 



2 

C**(R„R,) = d2fce^''-(^-^^)[l + T,(k)]a« (k) 



where 



T,(k)=tanhf:^^ 
' \2kBTJ 

Let us use the notation 

C^\p = C^^(Rj,Ri) = C$$(Rj,Ri±o) 

By using the definition (2.23) for the spectral function we easily obtain the following expressions 



Qr, = (^(i)»?Hi)) = Bo 
h 



1 — 77 

1 + Gq-UFq + ——Bo 



hi: 



U2 



r — — 



1 + Go + UFo 



1 — n 



Bo 



where we have defined 



= (r(^)€n*)) = I 

— ~^ Gi + UFi 

— 

Gn 
I > = 



Gi - UF, + Vt^Bi 



1 -n 



-Bi 



[a(k)]" 



2(27r)2 

n 

2(2^ 



/ d2^[a(k)]"[Ti(k)+r2(k)] 
JQb 



In the text, in order to compare different methods, we are also using the following correlation functions 



Ccc = (c(^)ct(^)) = + 2Q^ + 

= (c"(i)ct(z)) = + 2C-^ + C,, 



(c(z)77t(i)) = Q^ + C, 
Ce% = (c"(i)ryt(i)) = Cj« +C,, 



(A.2) 



(A.3) 



(A.4) 



(A.5) 



(A.6) 



(A.7) 



(A.E 



APPENDIX B: CONSERVATION OF THE FIRST FOUR SPECTRAL MOMENTA 



Let us show that the first four momenta of the electron spectral density satisfy the exact relations given by Eq. 
(2.26). By straightforward calculations we can see that: 
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M(0)(k)=F.T.({c(z),ct(j)}) = l 

M«(k) =F.T.(|(^z|) c(i),ct(j)|) = -M+fc/-4to(k) 



(B.l) 
(B.2) 



dt 



c{i),c\j) i> = -MM(i)(k) + f/ [m(k) - - [/)- - 4t(A + a(k)p) 



M(=^)(k) =F.T.( 

-4tF.T.J2M{{^§l)c{l),c\j)'^) = 
= At^ - yunC/ +^U^ + 16i^a^(k) + 8t/ia(k) - 8t|f/Q!(k) 
M(3)(k) = i^.T.(| (-^^7) 

On the other hand, from Eqs. 2.23 - 2.25 we have 

M(°)(k) = ^ [a«(k) + 2ag(k) + aW(k) 

1=1 

M(i)(k) = [a«(k) + 24; (k) + aW(k)] iJ,(k) 

'i=i 

= i?(k)M(°)(k) + g(k) 



dt 



(B.3) 



[-4ia(k) - fxf + 3-U [-4ta(k) - + -C/^ - 3-^[/2 - 8t-U^a{k) - AtAU^ - 4tpU^a{k) (B.4) 



= 1 



Ji-72 5(k) ^2.^(k) 



2 Q(k) Q(k) 



Prom Eq. (B.6) follows the identity 



h-h g(k) m(k) _ i?(k) M(i)(k) 

2 Q(k) + Q(k) g(k) ^ g(k) 



that leads to 



M(2)(k) = ^ [a«(k) + 2ag(k) + a«(k)] Ef{k} 

i=l 

= [R\k} + Q2(k)] MW(k) + 2i?(k)Q(k) 



i?(k) Af(i)(k) 



0(k) 0(k) 



At^ - fmU + ^-^U^ + 16t^a^{k) + 8tiia{k) - 8t^Ua{k) 



(B.5) 



(B.6) 



(B.7) 



(B.8) 



M(3)(k) = ^ [a«(k) + 24; (k) + a«(k)] £;f(k) 

= [i?3(k) +3i?(k)g2(k)] M(")(k) + [Q'{k)+3R\k)Q{k)] 
= 2ii(k)M(2)(k) + [Q2(k) - i?2(k)] M(i)(k) 



i?(k) ^ M(i)(k) 



g(k) Q(k) 



= [-4ia(k) - ij,f + [-4ia(k) - fj,f + -^U^ - Z-^^U^ - 8t jU^a{k) - AtMJ^ - AtpU^a{k) (B.9) 

Then, there is a complete equivalence among all considered theoretical approaches, when two-pole approximation is 
used. 
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